Hydrophobic interactions: an overview 
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We present an overview of tlie recent progress that has been made in understanding the origin 
of hydrophobic interactions. We discuss the different character of the solvation behavior of apolar 
solutes at small and large length scales. We emphasize that the crossover in the solvation behavior 
arises from a collective effect, which means that implicit solvent models should be used with care. 
We then discuss a recently developed explicit solvent model, in which the solvent is not described 
at the atomic level, but rather at the level of a density field. The model is based upon a lattice-gas 
model, which describes density fluctuations in the solvent at large length scales, and a Gaussian 
model, which describes density fluctuations at smaller length scales. By integrating out the small 
length scale field, a Ifamiltonian is obtained, which is a function of the binary, large-length scale field 
only. This makes it possible to simulate much larger systems than hitherto possible as demonstrated 
by the application of the model to the collapse of an ideal hydrophobic polymer. The results show 
that the collapse is dominated by the dynamics of the solvent, in particular the formation of a vapor 
bubble of critical size. Implications of these findings to the understanding of pressure denaturation 
of proteins are discussed. 

PACS numbers: 61.20.-p, 61.20.Gy, 68.08.-p, 82.70.Uv, 87.15.Aa 
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I. INTRODUCTION 

Hydrophobic interactions are widely believed to play a 
dominant role in the formation of large biological struc- 
tures 1^,^. Yet, the mechanism of the hydrophobic effect 
is still under debate. The solvation of small apolar species 
is well understood [^-^ . However, the attraction between 
two such species in water is weak [H^,^ and probably 
not responsible for the stability of biological structures. 
On the other hand, strong and long-ranged attractions 
have been measured between extended hydrophobic sur- 
faces |7|J^. But here, the origin of the effect is still be- 
ing discussed. It has been suggested that the interaction 
arises from electrostatic fluctuations , changes in water 
structure (lO|, bridging (sub)microscopic bubbles [0jl2) 
and a "drying" transition induced by the hydrophobic 
surfaces p^-[l5|. 

Recently, Lum, Chandler and Weeks have devel- 
oped a theory of hydrophobicity, which suggests that the 
nature of the hydrophobic effect precisely arises from 
the interplay of density fluctuations at both small and 
large length scales. Small apolar species only affect den- 
sity fluctuations in water at small length scales. Con- 
comitantly, water can only induce a relatively weak and 
short-ranged attraction between small hydrophobic ob- 
jects. In contrast, large hydrophobic species can affect 
density fluctuations at large length scales. At ambient 
conditions, water is close to phase coexistence. A suf- 
ficiently large hydrophobic object, or more importantly, 
an assembly of several small apolar species, can there- 
fore induce a depletion of water relative to the bulk den- 
sity . Recently, it has been demonstrated that this 
drying transition can induce a strong attraction between 
hydrophobic objects and provide a strong driving force 
for protein folding 



In this paper, we give an overview of the recent 
progress that has been made in understanding the ori- 
gin of the hydrophobic effect. In section |l| we discuss 
the statistics of density fluctuations at small and large 
length scales. Understanding these fluctuations is impor- 
tant, because it does not only provide insight into how 
the solvation free energy of solutes scales with their size, 
but also which models are needed to describe their sol- 
vation behavior. In the next section, we briefly discuss a 
recently developed model. In section IV, we see that this 
model gives a reasonable prediction of the solvation free 
energy of apolar solutes. In particular, it predicts that 
in the small length scale regime the solvation free energy 
scales with the size of the excluded volume of the solute, 
whereas in the large length scale regime it scales with 
the area of the excluded volume. The model also shows 
that the crossover in the solvation free energy arises from 
a collective effect in the solvent. Implicit solvent models 
cannot conveniently describe this collective effect. It thus 
appears that explicit solvent models are needed. Explicit 
atomistic solvent models, however, are computationally 
demanding. The model discussed in section III, lays the 



foundation for a scheme in which the solvent is not de- 
scribed at the atomic level, but rather at the level of a 
coarse-grained density field. This scheme, which is dis- 
cussed in section allows us to simulate the solvent 
much more efficiently. 

In section ^ we discuss the role of attractive inter- 
actions between the solutes and the solvent. Most of 
the theoretical work on the hydrophobic effect has fo- 
cused on the solvation behavior of ideal hydrophobic so- 
lutes |l^Jl^,^. These are objects that exclude solvent 
from a certain region in space, but have no attractive 
interactions with the solvent. In section VI, however, it 



is seen that the presence of weak dispersive interactions 
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does not significantly affect tfie solvation behavior of hy- 
drophobic solutes. 



Finally, in section VII we apply the scheme discussed 
in section ^ to study the collapse of an ideal hydropho- 
bic polymer. The simulations reveal that the dynamics 
of the collapse transition is dominated by the dynamics 
of the solvent. In particular, the rate limiting step is 
the formation of a vapor bubble of critical size. In ad- 
dition, we show that during the collapse the chain and 
the solvent remain out of equilibrium. Both observations 
imply that implicit solvent models should be used with 
great care. In the past, a statistical meaningful study of 
protein folding using explicit solvent models seemed im- 
practical. The current analysis, however, suggests that 
such studies become feasible by using a statistical field 
model like the one presented in sections 01 and 0. 



II. DENSITY FLUCTUATIONS AT SMALL AND 
LARGE LENGTH SCALES 

A theory of hydrophobicity should be able to correctly 
predict the solvation free energy of apolar solutes. The 
solvation free energy of a hard sphere, an example of an 
ideal hydrophobic object, is given by 



A/i(fcx) = -kBTln[P{N = 0;wcx)]. 



(1) 



Here ksT is Boltzmann's constant times temperature 
and P{N; Vcx) is the probability of finding N solvent 
particles inside a region, of volume Uex, from which 
the solvent would be excluded by the solute; note that 
P{N = 0;?;ex) corresponds to the insertion probability 
of a hard sphere. Computer simulations have revealed 
that density fluctuations that are entropic in origin, obey 
Gaussian statistics. For the hard-sphere fluid, the dis- 
tribution P{N; Vcx) (as a function of N) is found to be 
almost exactly Gaussian for spherical volumes Wcx that 
have a diameter d up to four times the diameter of the sol- 
vent molecules, a (i.e. hard-sphere diameter) [^. Con- 
ceivably, the distribution functions remain Gaussian for 
volumes that are larger than those that could be stud- 
ied in the simulations, albeit even larger voids could in- 
duce the formation of a solid-fluid like interface if the 
fluid is close to freezing. On the other hand, for vol- 
umes in SPC/E water at ambient conditions |6| and a 
Lennard-Jones liquid close to phase coexistence with the 
vapor the distribution P{N;Vcx) remains Gaussian 
up to approximately d = 1.0 — 2.0a. For such small 
volumes, the occupation statistics is dominated by en- 
tropic effects. The solvation free energy of small solutes 
is determined by the entropic cost of permitting only 
those density fluctuations that do not violate the effect 
of the strong forces that are present in the system. In 
a simple fluid, these strong forces arise from the repul- 
sive interactions between the particles; in water, these 
are the hydrogen bond forces. As the occupation statis- 
tics is Gaussian for small volumes, the insertion proba- 
bility P{N = 0; Wcx), and hence the solvation free energy 



of small solutes, can be obtained from a Gaussian the- 
ory. Examples of such Gaussian theories are the Pratt- 
Chandler theory of hydrophobicity P,p| and the theory 
of Hummer, Pratt and coworkers [|6|,p3[T 

For larger volumes the situation changes significantly. 
In a simple fluid, the attractive forces between the solvent 
molecules no longer cancel each other when the volume 
Wcx is excavated. The unbalanced attractive forces give 
rise to an unbalancing potential, which tends to push the 
solvent molecules away from the void. For water, a large 
excavated volume will drastically disrupt the water struc- 
ture. Close proximity of water to the void is energetically 
unfavorable, because the hydrogen bond network can no 
longer be maintained close to the surface of the void, as 
pointed out by Stillinger in 1973 |jl^. In both cases, sol- 
vent density is depleted near the surface of the void. This 
drying transition is a collective effect and can be inter- 
preted as a microscopic manifestation of a phase transi- 
tion. A Gaussian model cannot support such a transition 
as it only describes density fluctuations at microscopic 
length scales. For larger solutes, however, coarse-grained 
models, such as a Landau-Ginzburg model or a lattice- 
gas model, become useful. Such coarse-grained models 
describe density fluctuations at length scales larger than 
the bulk correlation length and can support phase transi- 
tions and sustain gas- liquid interfaces. The essence of our 
model, discussed in section [II, is to combine a Gaussian 



model with a lattice-gas model. 



III. THEORY OF SOLVATION 

Water at ambient conditions is a cold liquid. A cold 
liquid is a fluid that is well below the critical temperature. 
For such a cold liquid, the density of the vapor is typically 
orders of magnitude lower than that of the liquid. It is 
then natural to divide space into a (cubic) grid of cells 
and, taking the grid spacing, I, to be on the order of 
the bulk correlation length, ^, describe the states of the 
cells by a binary held, n.;, where = 1 if cell i is filled 
with "liquid" and = if cell i contains "vapor" . 
This binary field describes density fiuctuations at length 
scales larger than the grid spacing and it can support 
phase transitions. Density fiuctuations at length scales 
smaller than the grid spacing are described by a second 
field, 6p. This field does not support phase transitions, 
but it can describe the microscopic granularity of matter. 
As discussed above, we can assume that this field obeys 
Gaussian statistics. Hence, the full density p(ri) at point 
Ti is decomposed as: 



p(ri) = piHi + Sp{r^), 

where pi is the bulk liquid density. 

The partition function for our model is 



(2) 
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- = E / T^5p{v)C[{nk},5p{v)] 
X exp(-/3i/[K},5p(r)]), 



(3) 



where J 1)6 p{r) — J IliT>Sp{ri) denotes the functional in- 
tegration over the small length scale field, H[{nk} , 6p{r)] 
is the Hamiltonian as a functional of both Ui and Sp(r), 
and /3~^ is Boltzmann's constant times temperature, 
ksT. The quantity C[{nk} , 5p{r)] is a constraint func- 
tional. It has unit weight when the field 6p{r) together 
with {rii} satisfy whatever constraints are imposed by 
strong forces, and it is zero otherwise. Since {rii} and 
Sp(r) have greatly different character, the summation and 
integration in Eq. (^) do not redundantly count configu- 
ration space to any significant degree. 

The Hamiltonian of our model is given by the following 
expression: 



Hj - 1 



H[{nk} , Spir)] = HL[{nk}] - e' ^ J dr., Sp{r.,, ^^^^ 

i,j(nni) 



X 5p{r^) X \ri,Y'-{nk}] 5p{y'^) 

+ i?„o™[{nfc}]. (4) 

In the above expression, is the Hamiltonian of a 
lattice-gas model. It is given by 



HL[{nk}] 



(5) 



with p the imposed chemical potential, and e the inter- 
action parameter between nearest- neighbour cells. The 
lattice-gas model describes density fluctuations at length 
scales larger than the grid spacing and it can support 
phase transitions. Importantly, the energetic cost of cre- 
ating an interface is determined by the value of e. 

The term in Eq. ^ that is quadratic in Sp ensures the 
Gaussian weight of the small length scale field. In par- 
ticular, the response function x^^l^ii r^; {fifc}] determines 
the variance of this field. It depends upon the state of 
the large length scale field. If = 1 for all i, then 
X~^[ri, r^; {rifc}] reduces to the response function of the 

bulk liquid x~^[ri, Tj-; p;]- In general, we employ the ap- 
proximation that Sp{ri) = whenever rii = 0. Thus, 
Sp{r) is a Gaussian field, with a weight functional be- 
ing that of the bulk liquid, but constrained to be zero 
whenever rii = pc| ]. 

The second term on the right-hand-side of Eq. |4| is the 
term that couples the fluctuations of the field rii with 
those of the field Sp{r). As discussed by Lum et al. and 
Weeks and coworkers , for simple fluids the param- 

eter e' is related to the energy density of the bulk liquid. 
The last term on the right-hand-side of Eq. ^, -ffnorm, is a 
renormalization term, which ensures that, if the field 6p 



is integrated out in the absence of solutes, the Hamilto- 
nian iJ [{nfc}; 5p\ reduces to that of the lattice-gas model, 
HL[{nu]]. 

In the presence of a hard sphere, the constraint func- 
tional in Eq. is given by 



C[{nfc};(5p(r)] = [| 5[n,pi + Spiv,)]. 



(6) 



By exploiting the Fourier representation of delta func- 
tions, we can integrate out the field 5p{y). As discussed 
in detail in Ref. ||2^, we then arrive at the following ex- 
pression for the Hamiltonian: 



H[{7ik}]^HL[{nk}] 
+ kBT 

i.'j(occ) 



rii [piVi + fi] [piVj + fj] 



+ ksT In ■y/27r(T^; 
HL[{nk}]+Hs[v^^;{nk}\. 



(7) 



Here, 



= / drj{r,) ^n,v,€' Kj^ ^ (n^ - 1), (8) 



with K the isothermal compressibility, and 



c'^cx = / dr dr'x{r,r';pi). 



(9) 



The sum over j, j(occ) is over cells i and j that are occu- 
pied by the solute and Vi is the volume occupied by the 
solute in cell i [^ . 

The term Hs[vcx', {rik}] contains all the effects of the 
interaction between the ideal hydrophobic solute and the 
solvent. We stress that the interaction term solely arises 
from the constraint that is imposed upon the allowed 
density fiuctuations in the solvent. The excess chemical 
potential of the solute may be obtained by averaging this 
interaction free energy: 



f3Ap{v^ 



ln(exp(-/3i75[ucx; {nk}])L, 



(10) 



where denotes an ensemble average with the Hamil- 
tonian of the lattice gas, HL[{nk}] (see Eq. ||). 



IV. SCALING BEHAVIOR OF SOLVATION FREE 
ENERGIES 

Fig. 1^ shows the excess chemical potential of a hard 
sphere as a function of its radius, i?, in water at ambient 
conditions, as predicted by Eqs. ^ and For compar- 
ison, we also show the results of a Gaussian model and 
the results of a molecular simulation of a hard sphere in 
SPC/E water [||. 



3 



30.0 




— current model (Eq. 10) 
o-^ SPC/E 

-- Gaussian model (Eq. 12) 
■ ■ ■ Macroscopic model (Eq. 13) 



0.5 
R (nm) 

FIG. 1. Comparison of the results of the full model, based 
on Eq. ^ and Eq. |l^, with the predictions of the Gaussian 
model, Eq. and the results of a molecular simulation of a 
cavity in SPC/E water jl7| for the excess chemical potential 
per unit area of a hard sphere in water as a function of its size; 
the hard sphere excludes water from a spherical volume of ra- 
dius R. The horizontal line lies at the value of the surface ten- 
sion 7 of the vapor-liquid interface of water. The parameters 
of the full model are: I = 4.2A , e = 6.02fcsr, e' = 15.2fcsr, 
and n — picocx = 5.5 x 10~*fcsr. These parameters are cho- 
sen such that the surface-tension, the energy density and the 
imposed chemical potential of the lattice-gas model matches 
the corresponding values of water (see Ref. [po|). 

A. Small length scale regime 

In the small length scale regime, our model reduces 
to a Gaussian model. Small solutes only affect density 
fluctuations in the solvent at small length scales. As a 
consequence, rii — 1 for nearly all i. With Ui — 1 for all i, 
Hl [{iT'k}] and i?norm }] Eq. ^ become constants and 
thus irrelevant. Further, the response function reduces to 
that of the uniform fluid and the term that couples the 
density fluctuations of the small and large length scale 
field, becomes zero. As such, the Hamiltonian reduces to 
that of a Gaussian model ||] : 



Hg [Spiv)] 



dr I dr'Sp{r)x ^{r,r'; pi)6p{r'), 

(U) 



with 6p{r) = p{r) — pi, and the response function 
X^^(r, r';pi) being the response function of the uniform 
fluid. The Pratt-Chandler theory of hydrophobicity is 
based upon such a Gaussian model |^ . Similarly, we can 
directly obtain the solvation free energy of a Gaussian 
model from Eq. UU. It is given by: 



/3A/x(uox) ^ Pi Wcx/2o-t.cx + In V^TTov 



(12) 



Fig. 1^ shows that the full model and (thus) the Gaus- 
sian model give an excellent prediction for the solvation 
free energy of small solutes. The reason, as discussed in 
more detail in section ||, is that at small length scales, 
density fluctuations are highly Gaussian. A second point 
to note is that for small solutes, the excess chemical po- 
tential scales with the excluded volume of the solute. In 
this regime, the excess chemical potential is dominated 
by the entropic cost of constraining small length scale 
fluctuations, which is proportional to the volume of the 
space over which the solvent density fluctuations are con- 
strained. 



B. Large length scale regime 

For solutes larger than R > 4A, the prediction of the 
Gaussian model diverges from both the full model and 
the molecular simulation results. The divergence is due 
to a cavitation transition in the solvent. This cavitation 
transition can be understood as a microscopic liquid-gas 
phase transition in the solvent that is induced by the 
solute. In order to describe the solvation behavior of large 
hydrophobic objects in this "drying" regime, a theory of 
solvation should be able to support a liquid-gas phase 
transition and capture the enthalpic effects of creating 
a gas-liquid interface. The full model discussed in the 
previous section is such a model. A Gaussian model, 
on the other hand, cannot support a gas-liquid phase 
transition, nor sustain a gas-liquid interface, as it is based 
upon a density expansion around the uniform fluid. 

Fig. |l| shows that in the large length scale regime, the 
molecular simulation results for the excess potential are 
well described by 



(13) 



This is also the result for the solvation free energy of the 
theory of Hummer, Pratt and coworkers ||6|,p3[|. 



The first term on the right-hand side is the work to ex- 
pand a cavity against the external pressure, p. The sec- 
ond term on the right-hand side of the above equation 
describes the work to form an interface between the so- 
lute and the solvent, including a first-order "Tolman" 
correction due to the curvature of the interface. 

At ambient conditions, the pressure is very small, and 
the pV contribution to the free energy is only significant 
for extremely large solute volumes. Moreover, as dis- 
cussed in detail in Ref. |17|, 7 is well approximated by 
the surface tension of the vapor-liquid interface at coex- 
istence. Indeed, for large (but not unphysically large) 
solutes, the solvation free energy is dominated by the 
work to create a gas-liquid interface. 

It is often assumed that the solvation free energy of a 
hydrophobic species is proportional to the exposed sur- 
face area 26|. The molecular simulations, as well as 
the theoretical analyses, indicate that for solutes of bi- 
ologically relevant size, this is a reasonable assumption 
in the drying regime. For small solutes, however, the 
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excess chemical potential of a hydrophobic species does 
not scale with the area of the interface between the so- 
lute and the solvent. Rather, in the small length scale 
regime, the excess chemical potential scales with the ex- 
cluded volume of the solute. Finally, we point out that 
while macroscopic thermodynamic descriptions such as 
that underlying Eq. |l^ are useful for understanding the 
solvation behavior of large solutes (i? > 8 — 10 A), they 
are of limited use for describing the solvation behavior 
of small solutes. In the small length scale regime, micro- 
scopic models, such as the Gaussian model discussed in 
the previous section, are required. 



V. LANGEVIN DYNAMICS IN ISING SOLVENT 



The model discussed in section III can be used to de- 



velop a scheme in which only the solutes are treated at 
the atomic level; the solvent is described in terms of the 
binary, large-length scale density field Ui. In this way, 
the solvent can be simulated much more efficiently than 
using an explicit atomistic solvent model. The scheme is 
generally applicable to describe the motion of a collection 
of solutes, but here we will confine our attention to the 
motion of an ideal hydrophobic polymer consisting of Ng 
hard spheres. 

Before we discuss the scheme in detail, we point out 
that, in practice, it is reasonable to simplify the Hamil- 
tonian in Eq. ^ by neglecting the relatively small off- 
diagonal elements ||l8|. This yields the following Hamil- 
tonian for the coarse-grained density field: 

i?(r^=; {rifc}) = ^ (-/X + A/iox(wi))ni - e ^ niUj. 

i <i,j> 

(14) 



In the above expression, Vi = ^s,i, where Vs^i is 

the volume that monomer s occupies in cell i and Ng is 
indeed the number of monomers with coordinates . 
Further, A/iex(i'i) is the reversible work to accommodate 
a volume Vi in the solvent. We have taken A/icx to be 
proportional to the excluded volume: 



A^ox(Wi) « CVi 



(15) 



with c = 2.67 X lO^J/m'^ = GBkBT/nm"^ at room temper- 
ature. 

From H(r^'';{nk}), we can construct a free-energy 
functional, Q, for the field {uk} in the mean-field ap- 
proximation: 

f](r^=;K}) = 

J2 kBT{{n,) ln(n,;) + (1 - {n,)) ln(l - (n,))) 

i 

+ II (-A^ + Aa*cx(«»)) {n^) -\Y. (16) 

i <i,j> 



We can now construct a Car-Parinello scheme |2^,^3l to 
propagate the solvent. To this end, we define the La- 
grangian: 



i 

Here Vss (r^° ) is the intra-chain potential that describes 
the direction interaction between the monomers, M is 
the mass of the monomers and m is a fictitious mass 
associated with the dynamical variable {rii). From the 
Lagrangian, the following equation-of-motions for the 
monomers and the variables (n^) is obtained: 

Mr^ = -V^C = -Vo.[Vss{r^') + r!(r^^; {n,.})] (18) 



m{hi) = 



d{ni) 



d{ni) 



(19) 



Propagating the solvent field {rife } according to the above 
equation-of-motions with a proper choice of both the fic- 
titious mass m and the kinetic energy associated with 
{(rtfe)}, ensures that the free energy of the solvent is 
close to its minimum for each configuration of the poly- 
mer . This approach is similar in spirit to the method 
employed by Lowen, Madden and Hansen to simulate 
counterion screening in colloidal suspensions of polyelec- 
trolytes (29j. However, as we will demonstrate below, 
we cannot assume that the chain moves slowly on time 
scales for which disturbances in the solvent may relax. 
As we will see, the collapse of the polymer is driven by 
a cavitation transition in the solvent. Such a rare event 
cannot be captured by the above Car-Parinello scheme as 
nucleation barriers cannot be crossed. This means that, 
using a Car-Parinello scheme, the solvent would remain 
in its metastable state until it would become unstable. 
We therefore developed a novel scheme in which the dy- 
namics of the solvent and the chain are treated together. 

First, it should be realized that the small length scale 
field, has been integrated out. This means that tra- 
jectories are only true to nature on time scales larger than 
those required to relax the small length scale field ^p. 
This relaxation time is on the order of picoseconds , 
which implies that the dynamics is diffusive. We thus 
constructed a stochastic dynamics scheme. The elemen- 
tary step of the algorithm consists of propagating the 
polymer for Ms Langevin steps in the field of constant 
{n/c}, followed by a full Glauber sweep Q over the sol- 
vent field {rifc}. One Langevin step corresponds to: 

Ya{t + 5ts)= Va{t) 

+ — (- V„[ys5(r^'^ ) + (r^= ; {nu])] + <5F). (20) 
7 

where 7 is a friction coefficient, the value of which can be 
obtained from the diffusion constant of a single monomer 
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in water. The beads also experience a random force 5F; it 
is the dynamical remnant of the small length-scale field. 
In order to obtain a physically meaningful value for Mg, 
we have to compare the time step for the propagation of 
the polymer, Sts, to the time scale that corresponds to a 
Glauber sweep over the solvent variables. An estimate for 
the latter can be obtained by estimating the correlation 
time for a density fluctuation of length scale /: 6ti — 
l/[D(27r/Z)^], where D is the self-diffusion constant of 
liquid water. The value of Ms is thus given by Mg — 
6ti/6ts = 36. 

VI. ATTRACTIONS 

So far, we have discussed the solvation of ideal hy- 
drophobic solutes - objects that have no attractive inter- 
actions with the solvent molecules. Naturally, real hy- 
drophobic molecules have some affinity for water, and 
vice versa, because of the ubiquitous Van der Waals in- 
teractions. Huang and Chandler have studied in detail 
the effect of weak solute-solvent attractions on the solva- 
tion of non-polar molecules in water |^ . Their analysis 
was performed using an extension and improved parame- 
terization of the theory of Lum, Chandler and Weeks . 

Huang and Chandler showed that for small solutes, the 
presence of attractions has little effect on the solvent den- 
sity around the solutes, as has previously been appreci- 
ated theoretically |3| and as observed in simulations . 
The solvation behavior can be understood from the obser- 
vation that the attractive forces are very small in compar- 
ison with hydrogen bond forces, and their effects can be 
estimated by assuming that the water structure around 
a hard sphere is unaltered by adding an attraction to wa- 
ter. In particular, attractive interactions should produce 
a simple additive contribution to the solvation enthalpy, 
but no significant effect on the solvation entropy. 

For large hydrophobic species, attractions do have a 
notable effect. In the case of hard spheres, ie. without 
attractions, the solvent density near the surface of the 
solute is strongly depleted relative to that in the bulk 
liquid. In contrast, in the presence of attractions, the 
solvent density near the surface of the solute is close to 
that in the bulk liquid {g{R'^) « 1.24, where g{r) is the 
solvent radial distribution function and R is the radius 
of the hydrophobic sphere). This, however, is a result 
of the fact that the drying layer is very compressible. 
In particular, there is little free-energy cost to move the 
vapor-liquid interface. Due to this small energetic cost, 
the addition of an attractive potential as weak as that 
between alkanes and water is sufficient to draw the drying 
layer into contact with the hydrophobic surface. 

While drawn into contact with the hydrophobic sur- 
face, this interface is distinct from the interface that 
surrounds a small hydrophobic solute. The contact val- 
ues for g(r) for the large hydrophobic solutes are close 
to one, while the small hydrophobic solutes have con- 



tact values larger than two. More importantly, as for 
ideal hydrophobic objects, the solvation free energy scales 
with the size of the excluded volume for smaller solutes, 
whereas the solvation free energy scales with the area of 
the excluded volume for larger solutes. In addition, the 
crossover from the small to the large length scale regime 
is around 8 — lOA for both the ideal hydrophobic objects 
and the hydrophobic solutes that have attractive inter- 
actions with the solvent. Thus, the scaling behavior of 
the solvation free energy is not affected by the presence 
of weak attractions between the solute and the solvent. 
We therefore expect that weak dispersion forces will not 
have a strong effect on the role of hydrophobic forces in 
biological self-assembly. 

VII. HYDROPHOBIC POLYMER COLLAPSE 

Hydrophobic interactions have long been considered to 
play an important role in protein folding. We have per- 
formed computer simulations to study the collapse of a 
polymer consisting of Ns = 12 hard spheres in water. 
Even though this allows us to study the effect of hy- 
drophobic interactions on protein folding in perhaps its 
most basic form, the conventional approaches will not 
work. The theoretical analyses as discussed above, sug- 
gest that the hydrophobic effect arises from a collective 
effect in the solvent. Implicit solvent models cannot con- 
veniently capture this effect. It thus appears that ex- 
plicit solvent models are required. Explicit solvent mod- 
els, however, are computationally demanding. The anal- 
ysis is further complicated by the fact that the crossover 
from solvation in the small length scale regime to sol- 
vation in the large length scale regime occurs at around 
a nanometer. This implies that in order to study the 
hydrophobic effect, large system sizes are required. For 
the required system sizes, an explicit atomistic model of 
water would not be tractable in the simulations. In con- 
trast, the scheme discussed in section makes it possible 
to study the polymer collapse in water. 

Brute-force simulations, confirmed by the analysis dis- 
cussed below, suggest that the collapse transition is a rare 
event. We therefore performed transition-path-sampling 
simulations to harvest the rare, but representative, 
trajectories from the extended coil to the collapsed glob- 
ule state. Fig. || shows an example of such a trajectory 
at room temperature. Initially, the polymer is in the 
extended coil state. Then, by some spontaneous fiuctua- 
tion, the polymer partly collapses. The collapsed section 
forms a sufficiently large hydrophobic cluster that a va- 
por bubble is nucleated. Finally, the vapor bubble drives 
the beads of the polymer together. 
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A. Free energy of collapse 

Visual inspections of this and similar trajectories sug- 
gest that the collapse arises from an interplay between 
the size of the polymer and the formation of a bubble. 

We therefore mapped out the free-energy landscape as 
a function of the squared radius-of-gyration of the poly- 
mer, Rg, and the size of the largest vapor bubble present 
in the system, U. A contour plot of the free-energy 
landscape is shown in Fig. ^. It is seen that the path 
from the coil to the collapsed globule is one where, ini- 
tially, the radius-of-gyration decreases, while the size of 
the largest vapor bubble is still essentially zero. In this 
regime, the solvent still wets the polymer and the free 
energy hardly changes. When the radius-of-gyration be- 
comes small enough, however, a vapor bubble is formed. 
The free energy now sharply increases by some Ofc^T, 
until it reaches a saddle-point at {U,Rg) = (98,23.5Z^). 
From here on, the bubble grows spontaneously and drives 
all the beads of the polymer together. 

A close examination of the free-energy landscape re- 
veals that there is a small barrier of some 2kBT at 
{U,Rg) — (6.5, 70.0Z^), which separates the coil from 
a metastable intermediate. The presence of these two 
states arises from a competition between the entropy of 
the chain, which favors the fully extended coil state, and 
depletion forces, which favor the intermediate state. De- 
pletion forces are caused by the reduction in the volume 
from which the solutes exclude the solvent, when the so- 
lutes come together. The attraction between two small 
hydrophobic objects, which cannot induce a drying tran- 
sition in the solvent, predominantly arises from this ef- 
fect However, the free-energy landscape shows 
that this driving force is relatively small. The large driv- 
ing force for the collapse of the polymer comes from the 
drying transition. The free-energy difference between the 
fully-collapsed globule and the intermediate is SO/csT, 
whereas the free-energy difference between the interme- 
diate and the coil is only a few ksT. 

The collapse transition arises from an interplay be- 
tween density fluctuations of the solvent at small and 
large length scales. In the coil state, the monomers are 
well separated and the solvation free energy is dominated 
by the entropic cost of constraining small length scale 
fluctuations. This cost scales with the size of the ex- 
cluded volume and, to a good approximation, is given by 
^jA/Ltex(vi)- But when the monomers come together, 
this entropic cost is larger than the energetic cost of form- 
ing a vapor bubble that envelopes the monomers. As wa- 
ter is close to phase coexistence, this cost is dominated by 
the work to form an interface, which scales with the area 
of the excluded volume. Indeed, it is the crossover in the 
scaling behavior of the solvation free energy, as shown in 
Fig. 1^, which is the origin of the collapse transition. 

In order to make the above analysis more quantitative, 
we write the free energy of "folding" as: 

AGfold — AGsolv + AGintra- (21) 




FIG. 2. Four configurations from a trajectory where a 12 
unit hydrophobic chain in water goes from the coil to the glob- 
ule state, (a) shows a configuration from the equilibrated coil. 
The chain remained in configurations like that throughout a 
10 ns run at room temperature (T = 0.663e). On a much 
longer time scale, about 10~^ s, the chain typically does ex- 
hibit a transition from coil to globule. Such events are found 
with transition path sampling, equilibrating from an initial 
high temperature {T = 0.74e) 10 ns trajectory that exhibited 
the transition spontaneously. Three configurations from an 
equilibrated 1.5 ns trajectory that exhibits the collapse tran- 
sition at room temperature, are shown in (b), (c) and (d), 
with that in (c) being a configuration from the transition state 
surface. The transparent cubes denote the vapor cells. Those 
seen far from the chain are typical spontaneous density fiuctu- 
ations in bulk liquid water. The size of the simulation box is 
397 nm'^, corresponding to 42,875 cells. The potential energy 
function is given by Vss{r'^')+H{r^' ; {uk}). The latter term 
describes the solvent and the interaction between the polymer 
and the solvent; it is given by Eq. [14. The parameters are 
/ = 0.21 nm, ^1 - /icocx = 2.25 x lO'^sT, and e = 1.51fcsT. 
The intra-chain potential, Vss{r^"), is a function of the po- 
sitions of the centers of each of the 12 hydrophobic spheres 
(the red particles in the figure). It contains three parts: (1) 
steep (essentially hard sphere) repulsions between solute par- 
ticles such that their interparticle separations are larger than 
a — 0.72 nm; (2) stiff harmonic potentials bonding adjacent 
particles in the hydrophobic chain, ifcs((T— |ra+i — r^l)^, with 
ks = 14.1J/m^; (3) a bending potential favoring an extended 
chain, ^kgOa^, where 9a is the angle between (ra+2 — fa+i) 
and (rc+i — re), and ko = 1.85 x 10~^°J/rad^. The volumes 
Vi excluded from water by the chain are dynamic as they 
change with changing chain configuration, i.e., Vi = Vi ({r^}). 
Specifically, these volumes are computed by assuming water 
molecules have van der Waals radii equal to 0.14nm, and that 
the diameter of each hydrophobic unit is cr = 0.72 nm.. That 
is, points in the excluded volume, r, are those in the union of 
all volumes inscribed by |r — rc,| < 0.5 nm, a = 1, 2, Ns- 

Here AGsoiv is the reversible work to transfer the chain 
from its "wet" solvated state to the "dry" state in the va- 
por bubble, and AGintra denotes the free-energy change 
associated with the internal degrees of freedom of the 
polymer (see Fig. ^). In order to estimate AGsoiv, we 
make the following assumptions: (1) the fully collapsed 
globule is confined in a spherical bubble of radius R; (2) 
the volume, from which the polymer excludes solvent in 
the extended coil state is equal to the volume of the bub- 
ble that contains the globule; (3) in the extended coil 
state, the solvent wets the polymer and rii — 1 for all i. 
We then arrive at the following expression for AGsoiv: 
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AGsoiv = 47ri?2^ + -t:R^p - ^ A^la^{v^) (22) 



AttR^j + --kR^{p - c), 
o 



(23) 



where 7 is the surface tension of the vapor-liquid inter- 
face, p is the external pressure and c is given by Eq. 
Let us now compare the following numbers. First, from 
the measured free-energy difference between the coil and 
the globule, AGfoid ~ — SOfc^T, the radius of the bub- 
ble, R ~ 1.1 nm, and Eqs. |l| and |2^, we can obtain 
AGintra- This yiclds AGintra « TSfcsT. This should be 
compared with the respective contributions to AGsoiv, 



47ri?2^ « 258/cbT and 



SGlfefiT, which shows 



that the contribution of the internal degrees of freedom 
of the polymer to the free energy of collapse is rela- 
tively small. Indeed, the free energy of collapse of the 
hydrophobic polymer is dominated by the solvation free 
energy. Further, at ambient conditions p « c, which 
shows that the collapse transition is dominated by a com- 
petition between the energetic cost of creating a vapor- 
liquid interface in the globule and the entropic cost of 
maintaining the polymer in the wet state in the coil. 

This simple model also allows us to address the exper- 
imental observation that a wide range of proteins dena- 
ture under pressure. If we assume that (1) the specific 
volume of folding is determined primarily by AGsoiv and 
that (2) the dependence of the surface free energy upon 
pressure can be neglected, then we arrive at: 

A-Wfoid ~ — — ~ Vil - — ). (24) 

To our knowledge, no values for the dependence of the 
solvation free energy upon pressure for hard spheres of 
the size used here, have been reported. We therefore 
make some drastic assumptions. In order to compute 
dc/dp we assume that c « k/ kt, where fc is a constant, 
and K, is the isothermal compressibility. With ht ~ 4.6 x 
10-i°Pa"^ at 1 atmosphere and kt « 3.5 x lO-^^Pa"^ 
at 1000 atmosphere, we arrive at dc/dp w 0.82. This 
yields Awfow = Vil. - 0.82) « GOOcm^/mol. This is in 
agreement with the experimental observation that pro- 
teins denature under pressure, although values reported 
in the literature for proteins of similar size are signifi- 
cantly smaller, by about a factor six psIJS^. 



B. Dynamics of collapse 

The free-energy landscape shows that the rate-limiting 
step is the formation of a vapor bubble of critical size. 
To investigate the extent to which R^ and U also cor- 
rectly describe the dynamics of this nucleation event, 
we performed extensive transition-path-sampling simu- 
lations on an ensemble of trajectories, each of length 
150 ps [l|]. A scatter plot of the trajectories is shown 
in Fig. BT To identify the transition state surface, we 




FIG. 3. Contour plot of the free energy landscape for the 
collapse of our hydrophobic polymer, computed by Monte 
Carlo umbrella sampling ji^, using the weight functional 
exp [-PiVssir'^') + H{v'^=;{nk}))\ . The curves of constant 
free energy are drawn as a function of the squared radius of 
gyration of the polymer, Ft^, and the size of the largest bubble 
in the system, U . Here we have used the "cluster" criterion 
that two vapor cells that are nearest neighbours, belong to the 
same bubble; the size of the bubble is given by the number 
of vapor cells. Neighboring lines of the free-energy landscape 
are separated by 2.5fcsT. Superimposed is a scatter plot (in 
black) of the harvested 150 ps trajectories going from the coil 
to the globule state. The transition states are indicated in 
green. The harvesting was performed with transition path 
sampling, making 8,400 moves in trajectoryspace, of which 
75% were shooting and 25% were shifting We find that 
the plateau regime of the flux correlation function is reached 
after 50-70 ps po|j34| implying that the typical commitment 
time for trajectories to pass over the barrier is of the order of 
O.lns. Given this time and the fact that the figure shows the 
free energy barrier separating the extended coil and compact 
globule states to be about QksT, the half life of the extended 
chain is about 0.1 ns x exp (9) ~ 10"^ s. 



have identified the configurations on each trajectory from 
which newly initiated trajectories have equal probability 
of landing in either the coil or globule . These con- 

figurations are members of the transition-state-ensemble. 
We project them onto Fig. ^. First, it is seen that the 
ensemble is located close to the saddle-point in the free- 
energy landscape. It should also be noted that the transi- 
tion state ensemble is slightly tilted in the {U, R?g) plane. 
This shows that the larger the polymer, the larger the 
size of the critical vapor bubble that has to be nucle- 
ated. However, the scatter of the transition-state ensem- 
ble from a line in the (U,R^) plane is notable, which 
indicates that at least one other variable in addition to 
U and Rg plays a pertinent role in the collapse transition. 

It is also noteworthy that the transition paths do not 
follow the lowest free-energy path from the metastable 
intermediate to the fully collapsed globule (Fig. ||). In 
particular, the transition states are not centered around 
the saddle point and the dynamical paths do not fol- 
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FIG. 4. Sketch of the coUapse transition of the hydropho- 
bic polymer. The patterned region depicts the solvent, water. 
The white region indicates the volume from which the poly- 
mer excludes solvent. 



low the steepest descend path from the saddle point to 
the fully collapsed globule. The reason is that the poly- 
mer and the solvent move on different time scales. The 
polymer moves on time scales of nanoseconds, while the 
solvent moves on time scales of picoseconds. As a result, 
when a vapor bubble is nucleated, the chain cannot re- 
spond on the time scale at which the bubble is formed. 
Only when the bubble has reached a size of around 400, 
does the polymer fully collapse into the globule state. 

The above observation indicates that free-energy sur- 
faces should be interpreted with care. The dynamical 
pathways for a transition are not always fully determined 
by the underlying free-energy surface. Fig. ^ shows that 
other aspects can be important. In fact, when we artifi- 
cially force the polymer to move a hundred times faster 
with respect to the solvent (i.e. Ms = 3600 instead of 
36 - see section then the transition paths do follow 
the lowest free-energy path from coil to globule. This 
clearly demonstrates that the dynamical trajectories are 
not only determined by the underlying free-energy sur- 
face, but that the natural dynamics of the system, or the 
dynamics of the algorithm, can be important as well. It 
also means that the transition state ensemble does not 
need to coincide with the dividing surface in the free- 
energy landscape, which is the surface that yields the 
highest transmission coefficient 



VIII. CONCLUSIONS 



effects. In this regime, large apolar species can induce 
a cavitation transition in the solvent. More importantly, 
in the small length scale regime, the solvation free en- 
ergy scales with the volume of the solute, whereas in the 
drying regime, the solvation free energy scales with the 
exposed area of the solute. This behavior is not only ob- 
served for ideal hydrophobic objects, but also for solutes 
that have weak dispersive interactions with the solvent. 
The crossover behavior of the solvation free energy from 
the wetting regime to the drying regime is important, 
because it could be of significance to the formation of bi- 
ological structures. In most biological systems, the size 
of the hydrophobic species is such that individual species 
are in the wetting regime, while assemblies of such species 
are in the drying regime. In the wetting regime, water 
can only induce a relatively weak attraction between two 
small apolar species. When several of these species come 
together, however, water can induce a strong attraction 
between them. 

A clear example of this process is given by the col- 
lapse of the hydrophobic polymer. The strong driving 
force for the collapse of the polymer, is provided by the 
cavitation transition in the solvent. Importantly, the 
cavitation transition is only induced when a sufficiently 
large number of monomers comes together. This means 
that implicit solvent models, in which the interactions 
between the solutes is described by a sum of two body 
terms, cannot describe the cavitation transition. Equally 
important, the dynamics of the hydrophobic collapse is 
dominated by the dynamics of the solvent, especially the 
formation of a vapor bubble near the surface of a nucleat- 
ing cluster of hydrophobic monomers. As implicit solvent 
models cannot capture the dynamics of the solvent, they 
arc of limited use in studying the effect of hydrophobic 
interactions on the dynamics of biological self-assembly. 

The connection between hydrophobic collapse and the 
cavitation/drying transition provides a simple explana- 
tion for both cold denaturation and pressure denatura- 
tion of proteins. The range and the strength of the in- 
teractions between the hydrophobic objects in the dry- 
ing regime, is smaller when the solvent is moved away 
from phase coexistence. Thus, the lowering of temper- 
ature and the increase of pressure destabilize hydropho- 
bic collapse, because both actions move the solvent away 
from phase-coexistence. Finally, we believe that fluores- 
cence resonance energy transfer (FRET) experiments 1 3q] 
and nuclear magnetic resonance (NMR) spectroscopy 1 3^] 
should make it possible to observe the drying transition 
in protein folding experimentally. 



In this paper, we have given an overview of the sol- 
vation of hydrophobic solutes. Its character is very dif- 
ferent at small and large length scales. At small length 
scales, solvation is dominated by entropic effects and the 
solvent still wets the surface of the solute, even when 
the solute is highly hydrophobic In contrast, at 

large length scales, solvation is dominated by enthalpic 
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